#============================================================================================================================================================================================#
##Replication file for Banners, Barricades, and Bombs: The Tactical Choices of Social Movements and Public Opinion
rm(list=ls())
#---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
ipak <- function(pkg){  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
                        if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
                        sapply(pkg, require, character.only = TRUE)
}

packages <- c("memisc", "car") 
ipak(packages)

setwd("/Users/connorhuff/Dropbox/HuffKruszewska/Bargaining Tactics Experiment/Scripts/Final Replication Materials")
setwd("/Users/dominikakruszewska/Dropbox/HuffKruszewska/Bargaining Tactics Experiment/Scripts/Final Replication Materials")

##Loading the Data
data <- read.csv("BBB.csv")

###Loading the Data
#data <- as.data.set(spss.system.file('BBB.sav'))
#data <- as.data.frame(data)
#head(data)

#============================================================================================================================================================================================#
#===============================================================Prepping the data for analysis===============================================================================================#
#============================================================================================================================================================================================#
##renaming the variables
##Do you think the government should negotiate, answers are yes/no
colnames(data)[which(colnames(data)=="q18")] <- "should_negotiate"

##recoding yes to 1 and no to 0
data$should_negotiate_numeric <- recode(data$should_negotiate, "c('Yes')=1; c('No')=0",as.factor.result=FALSE) #recoding factor variable to number
data$should_negotiate_numeric <- as.numeric(data$should_negotiate_numeric)

##For individuals that think the government should negotiate, how strongly do they feel
colnames(data)[which(colnames(data)=="q19")] <- "yes_negotiate_strength"
##For individuals that think the government should not negotiate, how strongly do they feel
colnames(data)[which(colnames(data)=="q20")] <- "no_negotiate_strength"

##Once the government and organization enter negotations, what should be the outcome of these (regional autonomy, etc)?
colnames(data)[which(colnames(data)=="q22")] <- "negotiated_outcome"
data$negotiated_outcome_numeric <- recode(data$negotiated_outcome, "c('Independence for the region. The region will establish a new democratic government that will be responsible for overseei')=3; c('Regional autonomy. The regional government will control its own economic affairs but the region will remain part of the')=2;
                                          c('No regional autonomy. The central government will maintain full control over the region')=1",as.factor.result=FALSE)
data$negotiated_outcome_numeric <- as.numeric(data$negotiated_outcome_numeric)

##constructing the four point scale for should the government negotiate
data$yes_negotiate_strength_numeric <- recode(data$yes_negotiate_strength, "c('Very strongly')=4; c('Not very strongly')=3",as.factor.result=FALSE)
data$no_negotiate_strength_numeric <- recode(data$no_negotiate_strength, "c('Very strongly')=1; c('Not very strongly')=2",as.factor.result=FALSE)
data$no_negotiate_strength_numeric[is.na(data$no_negotiate_strength_numeric)] <- 0 ##Recoding NAs as 0 so that can collapse the two vectors into one
data$yes_negotiate_strength_numeric[is.na(data$yes_negotiate_strength_numeric)] <- 0
data$negotiate_strength_all <- data$yes_negotiate_strength_numeric + data$no_negotiate_strength_numeric #can add together since every variable either has its true value or is zero

##tax outcome question
colnames(data)[which(colnames(data)=="q23b_1")] <- "tax_outcome"


##The variable named test gives the different treatment categories
demonstration <- grepl("Dem", data$test)
data$treat_dem <- as.numeric(demonstration)

occupy <- grepl("Occ", data$test)
data$treat_occ <- as.numeric(occupy)

bombing <- grepl("Bom", data$test)
data$treat_bomb <- as.numeric(bombing)

foreign <- grepl("F", data$test)
data$treat_foreign <- as.numeric(foreign)

separatist <- grepl("Sep", data$test)
data$treat_sep <- as.numeric(separatist)

##creating two different datasets for the two vignettes
data.sep <- subset(data, data$treat_sep==1)
data.sm <- subset(data, data$treat_sep==0)



#============================================================================================================================================================================================#
#===============================================================Figure 1=====================================================================================================================#
#============================================================================================================================================================================================#
##Figure 1: Differince in proportion as specified in the pre-analysis plan for 

##Separatist movements
bomb.ppl.sep <- subset(data.sep, data.sep$treat_bomb==1)
dem.ppl.sep <- subset(data.sep, data.sep$treat_dem==1)
occ.ppl.sep <- subset(data.sep, data.sep$treat_occ==1)

diff.in.prop.sep <- prop.test(x=c(sum(bomb.ppl.sep$should_negotiate_numeric), sum(occ.ppl.sep$should_negotiate_numeric)), n=c(nrow(bomb.ppl.sep), nrow(occ.ppl.sep)))
diff.in.prop.sep #difference in proportions


##Social movements
bomb.ppl.sm <- subset(data.sm, data.sm$treat_bomb==1)
dem.ppl.sm <- subset(data.sm, data.sm$treat_dem==1)
occ.ppl.sm <- subset(data.sm, data.sm$treat_occ==1)

diff.in.prop.sm <- prop.test(x=c(sum(bomb.ppl.sm$should_negotiate_numeric), sum(occ.ppl.sm$should_negotiate_numeric)), n=c(nrow(bomb.ppl.sm), nrow(occ.ppl.sm)))
diff.in.prop.sm #difference in proportions

##differences in proportions comparing bombing and occupations for separatist organizations and social movements respectively
effect <- diff.in.prop.sep$estimate[1]-diff.in.prop.sep$estimate[2]
effect2 <- diff.in.prop.sm$estimate[1]-diff.in.prop.sm$estimate[2]


diff.in.prop.occ.dem.sep <- prop.test(x=c(sum(occ.ppl.sep$should_negotiate_numeric), sum(dem.ppl.sep$should_negotiate_numeric)), n=c(nrow(occ.ppl.sep), nrow(dem.ppl.sep)))
diff.in.prop.occ.dem.sep

diff.in.prop.occ.dem.sm <- prop.test(x=c(sum(occ.ppl.sm$should_negotiate_numeric), sum(dem.ppl.sm$should_negotiate_numeric)), n=c(nrow(occ.ppl.sm), nrow(dem.ppl.sm)))
diff.in.prop.occ.dem.sm

##differences in proportions comparing demonstrations and occupations for separatist organizations and social movements respectively
effect3 <- diff.in.prop.occ.dem.sep$estimate[1]-diff.in.prop.occ.dem.sep$estimate[2]
effect4 <- diff.in.prop.occ.dem.sm$estimate[1]-diff.in.prop.occ.dem.sm$estimate[2]


pdf("Figure_1.pdf")
par(mfrow=c(2,1))
par(oma=c(0,0,0,0))
plot(effect3,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(diff.in.prop.occ.dem.sep$conf.int[1],3,diff.in.prop.occ.dem.sep$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Occupation with Demonstration as Baseline",side=3,cex=1.4,line=1)
points(effect4,1,cex=3,pch=20)
segments(diff.in.prop.occ.dem.sm$conf.int[1],1,diff.in.prop.occ.dem.sm$conf.int[2],1,lwd=4)
axis(2,at=c(1,3),labels=c("SM","Sep"),cex.axis=1.2)
axis(1,at=c(seq(-0.4,0.4,by=0.1)),labels=c("-40","-30", "-20", "-10", "0", "10", "20", "30", "40"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)

plot(effect,3,cex=3,pch=20,xlim=c(-0.4,0.4),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments(diff.in.prop.sep$conf.int[1],3,diff.in.prop.sep$conf.int[2],3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Bombing with Occupation as Baseline",side=3,cex=1.4,line=1)
points(effect2,1,cex=3,pch=20)
segments(diff.in.prop.sm$conf.int[1],1,diff.in.prop.sm$conf.int[2],1,lwd=4)
axis(2,at=c(1,3),labels=c("SM","Sep"),cex.axis=1.2)
axis(1,at=c(seq(-0.4,0.4,by=0.1)),labels=c("-40","-30", "-20", "-10", "0", "10", "20", "30", "40"))
mtext("Percentage Point Change",side=1,cex=1.1,line=2.3)
dev.off()


#============================================================================================================================================================================================#
#===============================================================Figure 2=====================================================================================================================#
#============================================================================================================================================================================================#
##Figure 2: Difference in means with four point scale

###want to compare against the less extreme category
data.bomb.occ.sep <- subset(data.sep, data.sep$treat_dem!=1)

##regression comparing bombing to occupation for the separatist organization
reg.bomb.occ.sep <- lm(negotiate_strength_all ~ treat_bomb, data=data.bomb.occ.sep)
summary(reg.bomb.occ.sep) #signficant at the 0.1 level

data.bomb.occ.sm <- subset(data.sm, data.sm$treat_dem!=1)

##regression comparing bombing to occupation for the social movement
reg.bomb.occ.sm <- lm(negotiate_strength_all ~ treat_bomb, data=data.bomb.occ.sm)
summary(reg.bomb.occ.sm) #significant at the 0.05 level

#comparing bombing to occupation for separatist movements on a four point scale
effect5<-summary(reg.bomb.occ.sep)$coefficients[,1][2]
se5<-summary(reg.bomb.occ.sep)$coefficients[,2][2]

#comparing bombing to occupation for social movements on a four point scale
effect6<-summary(reg.bomb.occ.sm)$coefficients[,1][2]
se6<-summary(reg.bomb.occ.sm)$coefficients[,2][2]



##regression comparing occupation to demonstration
###want to compare against the less extreme category. 
data.occ.dem.sep <- subset(data.sep, data.sep$treat_bomb!=1)

##regression comparing occupation to demonstration for the separatist organization
reg.occ.dem.sep <- lm(negotiate_strength_all ~ treat_occ, data=data.occ.dem.sep)
summary(reg.occ.dem.sep) #no effect

data.occ.dem.sm <- subset(data.sm, data.sm$treat_bomb!=1)

##regression comparing occupation to demonstration for the social movement
reg.occ.dem.sm <- lm(negotiate_strength_all ~ treat_occ, data=data.occ.dem.sm)
summary(reg.occ.dem.sm) #no effect 

#comparing demonstration to occupation for separatist movements on a four point scale
effect7<-summary(reg.occ.dem.sep)$coefficients[,1][2]
se7<-summary(reg.occ.dem.sep)$coefficients[,2][2]

#comparing demonstration to occupation for social movements on a four point scale
effect8<-summary(reg.occ.dem.sm)$coefficients[,1][2]
se8<-summary(reg.occ.dem.sm)$coefficients[,2][2]


pdf("Figure_2.pdf")
par(mfrow=c(2,1))
par(oma=c(0,0,0,0))
plot(effect7,3,cex=3,pch=20,xlim=c(-0.8,0.4),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments((effect7-1.96*se7),3,(effect7+1.96*se7),3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Occupation with Demonstration as Baseline",side=3,cex=1.4,line=1.5)
points(effect8,1,cex=3,pch=20)
segments((effect8-1.96*se8),1,(effect8+1.96*se8),1,lwd=4)
axis(2,at=c(1,3),labels=c("SM","Sep"),cex.axis=1.2)
axis(1,at=c(seq(-0.8,0.4,by=0.2)),labels=c("-0.8", "-0.6", "-0.4", "-0.2", "0.0", "0.2", "0.4"))
mtext("Estimated Marginal Effect",side=1,cex=1.1,line=2.3)


plot(effect5,3,cex=3,pch=20,xlim=c(-0.8,0.4),yaxt="n",ylab="",ylim=c(0,4),xlab="",xaxt="n")
segments((effect5-1.96*se5),3,(effect5+1.96*se5),3,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Bombing with Occupation as Baseline",side=3,cex=1.4,line=1.5)
points(effect6,1,cex=3,pch=20)
segments((effect6-1.96*se6),1,(effect6+1.96*se6),1,lwd=4)
axis(2,at=c(1,3),labels=c("SM","Sep"),cex.axis=1.2)
axis(1,at=c(seq(-0.8,0.4,by=0.2)),labels=c("-0.8", "-0.6", "-0.4", "-0.2", "0.0", "0.2", "0.4"))
mtext("Estimated Marginal Effect",side=1,cex=1.1,line=2.3)
dev.off()



#============================================================================================================================================================================================#
#===============================================================Figure 3=====================================================================================================================#
#============================================================================================================================================================================================#
#Figure 3: What should they get during negotiations? (3 point scale)
##No regional autonomy is a 1; Independence is a 3

##regression comparing bombing and occupation
sep_reg_bomb_occ_negotiated_outcome <- lm(as.numeric(negotiated_outcome_numeric) ~ treat_bomb, data=data.bomb.occ.sep)
summary(sep_reg_bomb_occ_negotiated_outcome) #not significant
effect9<-summary(sep_reg_bomb_occ_negotiated_outcome)$coefficients[,1][2]
se9<-summary(sep_reg_bomb_occ_negotiated_outcome)$coefficients[,2][2]

##regression comparing occupation and demonstration
sep_reg_occ_dem_negotiated_outcome <- lm(as.numeric(negotiated_outcome_numeric) ~ treat_occ, data=data.occ.dem.sep)
summary(sep_reg_occ_dem_negotiated_outcome)
effect10<-summary(sep_reg_occ_dem_negotiated_outcome)$coefficients[,1][2]
se10<-summary(sep_reg_occ_dem_negotiated_outcome)$coefficients[,2][2]


pdf("Figure_3.pdf")
par(mfrow=c(2,1))
par(oma=c(0,0,0,0))
plot(effect10,1,cex=3,pch=20,xlim=c(-0.25,0.25),yaxt="n",ylab="", xlab="")
segments((effect10-1.96*se10),1,(effect10+1.96*se10),1,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Occupation with Demonstration as Baseline",side=3,cex=1.4,line=1)
mtext("Estimated Marginal Effect",side=1,cex=1.1,line=2.3)

plot(effect9,1,cex=3,pch=20,xlim=c(-0.25,0.25),yaxt="n",ylab="", xlab="")
segments((effect9-1.96*se9),1,(effect9+1.96*se9),1,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Bombing with Occupation as Baseline",side=3,cex=1.4,line=1)
mtext("Estimated Marginal Effect",side=1,cex=1.1,line=2.3)
dev.off()


#============================================================================================================================================================================================#
#===============================================================Figure 4=====================================================================================================================#
#============================================================================================================================================================================================#
##Figure 4: What should they get during negotiations with continuous tax measure outcome?

#comparing bombing to occupation
sep_reg_bomb_occ_tax_outcome <- lm(tax_outcome ~ treat_bomb, data=data.bomb.occ.sep)
summary(sep_reg_bomb_occ_tax_outcome)
effect11<-summary(sep_reg_bomb_occ_tax_outcome)$coefficients[,1][2]
se11<-summary(sep_reg_bomb_occ_tax_outcome)$coefficients[,2][2]

#comparing occupation to demonstration
sep_reg_occ_dem_tax_outcome <- lm(tax_outcome ~ treat_occ, data=data.occ.dem.sep)
summary(sep_reg_occ_dem_tax_outcome)
effect12<-summary(sep_reg_occ_dem_tax_outcome)$coefficients[,1][2]
se12<-summary(sep_reg_occ_dem_tax_outcome)$coefficients[,2][2]


pdf("Figure_4.pdf")
par(mfrow=c(2,1))
par(oma=c(0,0,0,0))
plot(effect12,1,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="")
segments((effect12-1.96*se12),1,(effect12+1.96*se12),1,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Occupation with Demonstration as Baseline",side=3,cex=1.4,line=1)
mtext("Estimated Marginal Effect",side=1,cex=1.1,line=2.3)

plot(effect11,1,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="")
segments((effect11-1.96*se11),1,(effect11+1.96*se11),1,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Bombing with Occupation as Baseline",side=3,cex=1.4,line=1)
mtext("Estimated Marginal Effect",side=1,cex=1.1,line=2.3)
dev.off()


#============================================================================================================================================================================================#
#===============================================================Figure 5=====================================================================================================================#
#============================================================================================================================================================================================#
##demonstration as the baseline for the comparison with bombing
data.bomb.dem.sep <- subset(data.sep, data.sep$treat_occ!=1)

sep_reg_bomb_dem_tax_outcome <- lm(tax_outcome ~ treat_bomb, data=data.bomb.dem.sep)
summary(sep_reg_bomb_dem_tax_outcome) #signficant at the 0.05 level
effect13<-summary(sep_reg_bomb_dem_tax_outcome)$coefficients[,1][2]
se13<-summary(sep_reg_bomb_dem_tax_outcome)$coefficients[,2][2]


pdf("Figure_5.pdf", height=3.5)
par(mfrow=c(1,1))
#par(fin=c(1,2))
plot(effect13,1,cex=3,pch=20,xlim=c(-10,10),yaxt="n",ylab="", xlab="")
segments((effect13-1.96*se13),1,(effect13+1.96*se13),1,lwd=4)
abline(v=0,lty=3)
mtext("Marginal Effect of Bombing with Demonstration as Baseline",side=3,cex=1.4,line=1)
mtext("Estimated Marginal Effect",side=1,cex=1.1,line=2.3)
dev.off()






#============================================================================================================================================================================================#
#===============================================================Figure 6=====================================================================================================================#
#============================================================================================================================================================================================#
##Code for the Structural Topic Model
#---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
ipak <- function(pkg){  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
                        if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
                        sapply(pkg, require, character.only = TRUE)
}

packages <- c("stm", "tm", "wordcloud") 
ipak(packages)


#Loading the Data; Note: due to file encoding issues with other file formats and Polish need to load as a .dat file
data2 <- read.delim("BBB_Replication_Text.dat", fileEncoding="utf-8")
head(data2)


##Creating one column for all open ended responses
data2$all_text <- paste(data2$TQ21, data2$TQ23)


##Note: the variable named test gives the different treatment categories
##Create variables for treatments to use for stm analysis
data2$treat_sep_dem_fo <- ifelse(data2$TEST == 11, 1, 0)

data2$treat_sep_dem <- ifelse(data2$TEST == 12, 1, 0)

data2$treat_sep_occ_fo <- ifelse(data2$TEST == 13, 1, 0)

data2$treat_sep_occ <- ifelse(data2$TEST == 14, 1, 0)

data2$treat_sep_bom <- ifelse(data2$TEST == 15, 1, 0)

data2$treat_sep_bom_fo <- ifelse(data2$TEST == 16, 1, 0)

data2$treat_sm_dem <- ifelse(data2$TEST == 21, 1, 0)

data2$treat_sm_occ <- ifelse(data2$TEST == 22, 1, 0)

data2$treat_sm_bomb <- ifelse(data2$TEST == 23, 1, 0)

data2$bombing <- ifelse(data2$TEST == 15 | data2$TEST == 16 | data2$TEST == 23, 1, 0)


##processing the data - lower case, removes punctuation etc.

##stopwords: stopwords list from http://members.unine.ch/jacques.savoy/clef/
stoplist<-scan(file='/Users/connorhuff/Dropbox/HuffKruszewska/Bargaining Tactics Experiment/Data/polishST.txt', what='character', strip.white=TRUE, blank.lines.skip=TRUE)
stoplist<-scan(file='/Users/dominikakruszewska/Dropbox/HuffKruszewska/Bargaining Tactics Experiment/Data/polishST.txt', what='character', strip.white=TRUE, blank.lines.skip=TRUE)


stoplist <- c(stoplist, "sie", "aby", "żeby", "ale", "tym", "tego", "wiem", "zdania", "powinny", "powinni", "powinna", "powinien", "należy", "trzeba", "każdym", "będzie", "będziemy", "można", "uważam") #add frequently occuring, uninformative words

##words removed: "nie" (no)
stoplist <- stoplist[-83]


temp<-textProcessor(documents=data2$all_text,metadata=data2, stem=FALSE, removestopwords=FALSE, customstopwords=stoplist, language="na")


##preparing an object for stm analysis
out <- prepDocuments(temp$documents, temp$vocab, temp$meta)
documents <- out$documents
vocab <- out$vocab
meta <- out$meta


##choose three topics after examining output because after > 3, topics hard to distinguish

##model selection for k=3
mod.select <- selectModel(out$documents, out$vocab, K=3, prevalence=~bombing, data=meta, max.em.its=500, runs=50, seed=02138)

plotModels(mod.select) #should choose model on the upper right (based exclusivity and coherence)

bomModel <- mod.select$runout[[2]] #choose the second model

##examine output of the model, words associated with each topic to infer topic labels
summary.STM(bomModel)
labelTopics(bomModel, n=20)

##check representative documents for each topic as a validation step for labelling topics

#topic 1
thoughts1 <- findThoughts(bomModel,texts=meta$all_text, topics=1, n=1)$docs[[1]]
plotQuote(thoughts1)

#topic 2
thoughts2 <- findThoughts(bomModel,texts=meta$all_text, topics=2, n=1)$docs[[1]]
plotQuote(thoughts2)

#topic 3
thoughts3 <- findThoughts(bomModel,texts=meta$all_text, topics=3, n=1)$docs[[1]]
plotQuote(thoughts3)


##plot frequent and exclusive words for each topic in Polish
plot.STM(bomModel,type="labels", labeltype="frex", xlim=c(0,1), topic.names=c("Topic 1: Terrorism", "Topic 2: Negotiations and Agreement", "Topic 3: Citizen's Voice"))


pdf("Figure_6.pdf", width=8)
plot.STM(bomModel,type="labels", labeltype="frex", topic.names=c("Topic 1: Terrorism", "Topic 2: Negotiations and Agreement", "Topic 3: Citizen's Voice"), 
         custom.labels=c("negotiates, such, terrorists, know, movement, opinion, it's about, topic, separatists, others, let, talks, security, good, all, everyone, issue, terrorists, still, matters", 
                         "negotiations, agreement, difficult, tell, reach, agree, better, government, conversation, must, way, talks, sides, peace, country, problems, every, compromise, better, some", 
                         "law, people, nation, listen, saying, people, movement, situation, something, have, country, such, voice, lead to, citizens, goal, war, autonomy, democracy, right"))
dev.off()



#============================================================================================================================================================================================#
#===============================================================Figure 7=====================================================================================================================#
#============================================================================================================================================================================================#
##Figure 7: Effect of tactical on topic prevalence in open-ended responses

bomModel_effect <- estimateEffect(1:3 ~ bombing, bomModel, meta)


pdf("Figure_7.pdf")
plot.estimateEffect(bomModel_effect, covariate = "bombing", topics = 1:3, model=bomModel, method="difference", ci.level=0.95,
                    cov.value1="1",cov.value2="0", main="", xlab="Difference in Topic Proportions between Bombing and Baseline Categories",
                    xlim=c(-.2,.2), labeltype = "custom", custom.labels = c('Topic 1: Terrorism', 'Topic 2: Negotiations and Agreement', "Topic 3: Citizen's Voice"))

dev.off()



#============================================================================================================================================================================================#
#===============================================================Figure 8=====================================================================================================================#
#============================================================================================================================================================================================#
##Figure 8: Testing for Heterogeneous treatment effects for foreign

##Function used to estimate heterogeneous treatment effects for foreign
interaction_plot_binary2 <- function(model, effect, moderator, interaction, varcov="default", conf=.95, title="Marginal effects plot", xlabel="Value of moderator", 
                                     ylabel="Estimated Marginal Effect", factor_labels=c(0,1)){
  
  # Extract Variance Covariance matrix
  if (varcov == "default"){
    covMat = vcov(model)
  }else{
    covMat = varcov
  }
  
  # Extract the data frame of the model
  mod_frame = model.frame(model)
  
  # Get coefficients of variables
  beta_1 = model$coefficients[[effect]]
  beta_3 = model$coefficients[[interaction]]
  
  # Create list of moderator values at which marginal effect is evaluated
  x_2 <- c(0,1)
  
  # Compute marginal effects
  delta_1 = beta_1 + beta_3*x_2
  
  # Compute variances
  var_1 = covMat[effect,effect] + (x_2^2)*covMat[interaction, interaction] + 2*x_2*covMat[effect, interaction]
  
  # Standard errors
  se_1 = sqrt(var_1)
  
  # Upper and lower confidence bounds
  z_score = qnorm(1 - ((1 - conf)/2))
  upper_bound = delta_1 + z_score*se_1
  lower_bound = delta_1 - z_score*se_1
  
  # Determine the bounds of the graphing area
  max_y = max(upper_bound)
  min_y = min(lower_bound)
  
  # Initialize plotting window
  plot(y=c(), x=c(), xlim=c(min_y, max_y), ylim=c(-.5, 1.5), ylab=xlabel, xlab=ylabel, main=title, yaxt="n")
  
  # Plot points of estimated effects
  points(y=x_2, x=delta_1, pch=16)
  
  # Plot lines of confidence intervals
  lines(y=c(x_2[1], x_2[1]), x=c(upper_bound[1], lower_bound[1]), lty=1, lwd=4)
  lines(y=c(x_2[2], x_2[2]), x=c(upper_bound[2], lower_bound[2]), lty=1, lwd=4)

  # Label the axis
  axis(side=2, at=c(0,1), labels=factor_labels)
  
  # Add a dashed horizontal line for zero
  abline(v=0, lty=3)
  
}

##Note that the foreign variable only occurs when we have a separatist organization.
reg1 <- lm(should_negotiate_numeric ~ treat_bomb + treat_foreign + treat_bomb*treat_foreign, data=data.bomb.occ.sep)
summary(reg1)


pdf("Figure_8.pdf")
interaction_plot_binary2(reg1, effect="treat_bomb", moderator="treat_foreign", 
                         interaction="treat_bomb:treat_foreign", factor_labels=c("No","Yes"), 
                         xlabel="Incident in Foreign Country", 
                         title="Marginal Effect of Bombing by Location")
dev.off()









